#Replication materials for "Fed Up"
#Political Research Quarterly 2024
#Ian G. Anson
#iganson@umbc.edu
#2024-06-17

require(tidyverse)
require(ggthemes)
require(car)
require(stargazer)

#readin

dat <- read.csv("fed.csv")

glimpse(dat)

#clean scales from Qualtrics survey
dat$fed_rate <- car::recode(dat$fed_rate, "1=5;2=4;6=3;3=2;4=1")
dat$fed_power <- car::recode(dat$fed_power,"1=0;2=1;3=2;else=NA")
dat$fed_support <- car::recode(dat$fed_support,"1=5;2=4;3=3;4=2;5=1")
dat$rep1dem0 <- car::recode(dat$pid7, "1:3 = 0; 4 = NA; 5:7 = 1")

#rescale variables to run from 0:1 scales for easy interpretation
dat$fed_power2 <- (dat$fed_power)/2
dat$fed_rate2 <- (dat$fed_rate - 1)/4
dat$fed_support2 <- (dat$fed_support - 1)/4
dat$fed_learn <- car::recode(dat$fed_learn,"2=0")


#Fed performance evaluations

mod1 <- lm(data=dat,fed_rate2~as.factor(FL_9_DO))
summary(mod1)

mod2 <- lm(data=dat[dat$rep1dem0==1,],fed_rate2~as.factor(FL_9_DO))
summary(mod2)

mod3 <- lm(data=dat[dat$rep1dem0==0,],fed_rate2~as.factor(FL_9_DO))
summary(mod3)

stargazer(mod1,mod2,mod3,type="html",out="fedsupport.html",
          star.cutoffs=c(0.1,0.05,0.01))


#Perceptions that Fed is too powerful

mod4 <- lm(data=dat,fed_power2~as.factor(FL_9_DO))
summary(mod4)

mod5 <- lm(data=dat[dat$rep1dem0==1,],fed_power2~as.factor(FL_9_DO))
summary(mod5)

mod6 <- lm(data=dat[dat$rep1dem0==0,],fed_power2~as.factor(FL_9_DO))
summary(mod6)

#Support for Fed economic involvement

mod7 <- lm(data=dat,fed_support2~as.factor(FL_9_DO))
summary(mod7)

mod8 <- lm(data=dat[dat$rep1dem0==1,],fed_support2~as.factor(FL_9_DO))
summary(mod8)

mod9 <- lm(data=dat[dat$rep1dem0==0,],fed_support2~as.factor(FL_9_DO))
summary(mod9)

#Learning about the Fed

mod10 <- glm(data=dat,fed_learn~as.factor(FL_9_DO),
             family=binomial(link='logit'))
mod11 <- glm(data=dat[dat$rep1dem0==1,],fed_learn~as.factor(FL_9_DO),
            family=binomial(link='logit'))
mod12 <- glm(data=dat[dat$rep1dem0==0,],fed_learn~as.factor(FL_9_DO),
            family=binomial(link='logit'))

stargazer(mod10,mod11,mod12,type="html",out="fedlearn.html",
          star.cutoffs=c(0.1,0.05,0.01))



#manipulation checks

dat$fed_mcT2 <- grepl("1",dat$fed_mc)
dat$fed_mcT3 <- grepl("2",dat$fed_mc)
dat$fed_mcT4 <- grepl("3",dat$fed_mc)


modM1 <- lm(data=dat,fed_mcT2~as.factor(FL_9_DO))
summary(modM1)
modM2 <- lm(data=dat,fed_mcT3~as.factor(FL_9_DO))
summary(modM2)
modM3 <- lm(data=dat,fed_mcT4~as.factor(FL_9_DO))
summary(modM3)

stargazer(modM1,modM2,modM3, type="html",out="mc_checks.html",
          star.cutoffs=c(0.1,0.05,0.01))








###Graphics
#Fig. 1: Fed performance eval.

modelLabels <- c("Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment")

modelCoefs <- c(mod1$coefficients[2],
                mod1$coefficients[3],
                mod1$coefficients[4])

modelLow95 <- c(confint(mod1,level=0.95)[2,1],
                confint(mod1,level=0.95)[3,1],
                confint(mod1,level=0.95)[4,1]
                )

modelLow90 <- c(confint(mod1,level=0.90)[2,1],
                confint(mod1,level=0.90)[3,1],
                confint(mod1,level=0.90)[4,1]
)


modelHi90 <-  c(confint(mod1,level=0.90)[2,2],
                confint(mod1,level=0.90)[3,2],
                confint(mod1,level=0.90)[4,2]
)

modelHi95 <-  c(confint(mod1,level=0.95)[2,2],
                confint(mod1,level=0.95)[3,2],
                confint(mod1,level=0.95)[4,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$row <- c(1,2,3)
# 
# grafMat$pid <- c(rep("All Respondents",3),
#                  rep("Republicans Only",3),
#                  rep("Democrats Only",3))




graf1 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=modelLabels,x=modelCoefs),size=2)+
  geom_errorbarh(aes(y=modelLabels,xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1)+
  geom_errorbarh(aes(y=modelLabels,xmax=modelHi95,xmin=modelLow95),height=0.1)+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Experimental Treatments vs. Control, Fed Performance Rating")+
  xlab("Fed Performance Rating (0:1)")+
  ylab("")+
  #scale_x_continuous(limits=c(0,0.28),breaks=c(0,0.1,0.2))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=9))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  #scale_color_manual(values = c("black","#0072B2","#FF3333"))+
  theme(legend.position="none")

windows(8,6)
graf1

#Fig. 2: Fed too powerful?

modelLabels <- c("Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment")

modelCoefs <- c(mod4$coefficients[2],
                mod4$coefficients[3],
                mod4$coefficients[4])

modelLow95 <- c(confint(mod4,level=0.95)[2,1],
                confint(mod4,level=0.95)[3,1],
                confint(mod4,level=0.95)[4,1]
)

modelLow90 <- c(confint(mod4,level=0.90)[2,1],
                confint(mod4,level=0.90)[3,1],
                confint(mod4,level=0.90)[4,1]
)


modelHi90 <-  c(confint(mod4,level=0.90)[2,2],
                confint(mod4,level=0.90)[3,2],
                confint(mod4,level=0.90)[4,2]
)

modelHi95 <-  c(confint(mod4,level=0.95)[2,2],
                confint(mod4,level=0.95)[3,2],
                confint(mod4,level=0.95)[4,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$row <- c(1,2,3)

# grafMat$pid <- c(rep("All Respondents",3),
#                  rep("Republicans Only",3),
#                  rep("Democrats Only",3))




graf2 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=modelLabels,x=modelCoefs),size=2)+
  geom_errorbarh(aes(y=modelLabels, xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1)+
  geom_errorbarh(aes(y=modelLabels, xmax=modelHi95,xmin=modelLow95),height=0.1)+
  # facet_wrap(~modelLabels,ncol=1)+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Experimental Treatments vs. Control, Fed Too Powerful?")+
  xlab("Fed Too Powerful or Not Powerful Enough? (0:1)")+
  ylab("")+
  #scale_x_continuous(limits=c(0,0.28),breaks=c(0,0.1,0.2))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=9))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  # scale_color_manual(values = c("black","#0072B2","#FF3333"))+
  theme(legend.position="none")

windows(8,6)
graf2



#Fig. 3: Fed Legitimacy

modelLabels <- c("Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment",
                 "Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment",
                 "Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment")

modelCoefs <- c(mod7$coefficients[2],
                mod7$coefficients[3],
                mod7$coefficients[4],
                mod8$coefficients[2],
                mod8$coefficients[3],
                mod8$coefficients[4],
                mod9$coefficients[2],
                mod9$coefficients[3],
                mod9$coefficients[4])

modelLow95 <- c(confint(mod7,level=0.95)[2,1],
                confint(mod7,level=0.95)[3,1],
                confint(mod7,level=0.95)[4,1],
                confint(mod8,level=0.95)[2,1],
                confint(mod8,level=0.95)[3,1],
                confint(mod8,level=0.95)[4,1],
                confint(mod9,level=0.95)[2,1],
                confint(mod9,level=0.95)[3,1],
                confint(mod9,level=0.95)[4,1]
)

modelLow90 <- c(confint(mod7,level=0.90)[2,1],
                confint(mod7,level=0.90)[3,1],
                confint(mod7,level=0.90)[4,1],
                confint(mod8,level=0.90)[2,1],
                confint(mod8,level=0.90)[3,1],
                confint(mod8,level=0.90)[4,1],
                confint(mod9,level=0.90)[2,1],
                confint(mod9,level=0.90)[3,1],
                confint(mod9,level=0.90)[4,1]
)


modelHi90 <-  c(confint(mod7,level=0.90)[2,2],
                confint(mod7,level=0.90)[3,2],
                confint(mod7,level=0.90)[4,2],
                confint(mod8,level=0.90)[2,2],
                confint(mod8,level=0.90)[3,2],
                confint(mod8,level=0.90)[4,2],
                confint(mod9,level=0.90)[2,2],
                confint(mod9,level=0.90)[3,2],
                confint(mod9,level=0.90)[4,2]
)

modelHi95 <-  c(confint(mod7,level=0.95)[2,2],
                confint(mod7,level=0.95)[3,2],
                confint(mod7,level=0.95)[4,2],
                confint(mod8,level=0.95)[2,2],
                confint(mod8,level=0.95)[3,2],
                confint(mod8,level=0.95)[4,2],
                confint(mod9,level=0.95)[2,2],
                confint(mod9,level=0.95)[3,2],
                confint(mod9,level=0.95)[4,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$row <- c(1,1,1,2,2,2,3,3,3)

grafMat$pid <- c(rep("All Respondents",3),
                 rep("Republicans Only",3),
                 rep("Democrats Only",3))




graf3 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=factor(pid,level = c("Democrats Only","Republicans Only","All Respondents")),
                 x=modelCoefs, color=pid),size=2)+
  geom_errorbarh(aes(y=factor(pid,level = c("Democrats Only","Republicans Only","All Respondents")),
                     xmax=modelHi90,xmin=modelLow90, color=pid),height=0.0,size=1.1)+
  geom_errorbarh(aes(y=factor(pid,level = c("Democrats Only","Republicans Only","All Respondents")),
                     xmax=modelHi95,xmin=modelLow95, color=pid),height=0.1)+
  facet_wrap(~modelLabels,ncol=1)+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Experimental Treatments vs. Control, Perceptions of Fed Legitimacy")+
  xlab("Fed Legitimacy (0:1)")+
  ylab("")+
  #scale_x_continuous(limits=c(0,0.28),breaks=c(0,0.1,0.2))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=9))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values = c("black","#0072B2","#FF3333"))+
  theme(legend.position="none")


windows(8,6)
graf3


#desire for more info
#linearize
mod10g <- lm(data=dat,fed_learn~as.factor(FL_9_DO))

modelLabels <- c("Political Influence Treatment",
                 "Size and Scope Treatment",
                 "Agency Performance Treatment")

modelCoefs <- c(mod10g$coefficients[2],
                mod10g$coefficients[3],
                mod10g$coefficients[4])

modelLow95 <- c(confint(mod10g,level=0.95)[2,1],
                confint(mod10g,level=0.95)[3,1],
                confint(mod10g,level=0.95)[4,1]
)

modelLow90 <- c(confint(mod10g,level=0.90)[2,1],
                confint(mod10g,level=0.90)[3,1],
                confint(mod10g,level=0.90)[4,1]
)


modelHi90 <-  c(confint(mod10g,level=0.90)[2,2],
                confint(mod10g,level=0.90)[3,2],
                confint(mod10g,level=0.90)[4,2]
)

modelHi95 <-  c(confint(mod10g,level=0.95)[2,2],
                confint(mod10g,level=0.95)[3,2],
                confint(mod10g,level=0.95)[4,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$row <- c(1,2,3)


graf4 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=modelLabels,x=modelCoefs),size=2)+
  geom_errorbarh(aes(y=modelLabels,xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1)+
  geom_errorbarh(aes(y=modelLabels,xmax=modelHi95,xmin=modelLow95),height=0.1)+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Experimental Treatments vs. Control, Desire for More Info. About Fed")+
  xlab("Fed Info-Seeking (0:1)")+
  ylab("")+
  #scale_x_continuous(limits=c(0,0.28),breaks=c(0,0.1,0.2))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=9))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  #scale_color_manual(values = c("black","#0072B2","#FF3333"))+
  theme(legend.position="none")

windows(8,6)
graf4




##APPENDICES



#Descriptive Statistics

dat$male <- car::recode(dat$gender,"1=1;else=0")
dat$ideo <- car::recode(dat$ideo,"6=NA")

dat2 <- dat %>%
  select(age,male,race,edu,hhinc,employ,union,pid7,ideo,
         fed_rate,fed_power,fed_support,fed_learn,fed_mc,
         FL_9_DO)

stargazer(dat2, type="html", out = "descriptives.html")




#Subgroup analyses

#check 1: attentiveness

dat$check1 <- car::recode(dat$Q94, "'5,6'=1;else=0")
dat$check2 <- car::recode(dat$Q95, "3=1;else=0")
dat$check <- dat$check1 + dat$check2
dat$check <- car::recode(dat$check,"0:1=0;else=1")

mod1 <- lm(data=dat[dat$check==1,],fed_rate2~as.factor(FL_9_DO))
summary(mod1)

mod2 <- lm(data=dat[dat$check==1,],fed_power2~as.factor(FL_9_DO))
summary(mod2)

mod3 <- lm(data=dat[dat$check==1,],fed_support2~as.factor(FL_9_DO))
summary(mod3)


stargazer(mod1,mod2,mod3, type="html",out="robust1.html",
          star.cutoffs=c(0.1,0.05,0.01))


#political knowledge

dat$selfrate <- car::recode(dat$Q91,"1=5;2=4;3=3;4=2;5=1")

mod1 <- lm(data=dat[dat$selfrate>2,],fed_rate2~as.factor(FL_9_DO))
summary(mod1)

mod2 <- lm(data=dat[dat$selfrate>2,],fed_power2~as.factor(FL_9_DO))
summary(mod2)

mod3 <- lm(data=dat[dat$selfrate>2,],fed_support2~as.factor(FL_9_DO))
summary(mod3)


stargazer(mod1,mod2,mod3, type="html",out="robust2.html",
          star.cutoffs=c(0.1,0.05,0.01))


#check 3: ideology

mod1 <- lm(data=dat[dat$ideo>3,],fed_rate2~as.factor(FL_9_DO))
summary(mod1)

mod2 <- lm(data=dat[dat$ideo>3,],fed_power2~as.factor(FL_9_DO))
summary(mod2)

mod3 <- lm(data=dat[dat$ideo>3,],fed_support2~as.factor(FL_9_DO))
summary(mod3)


stargazer(mod1,mod2,mod3, type="html",out="robust3_libs.html",
          star.cutoffs=c(0.1,0.05,0.01))


mod1 <- lm(data=dat[dat$ideo<3,],fed_rate2~as.factor(FL_9_DO))
summary(mod1)

mod2 <- lm(data=dat[dat$ideo<3,],fed_power2~as.factor(FL_9_DO))
summary(mod2)

mod3 <- lm(data=dat[dat$ideo<3,],fed_support2~as.factor(FL_9_DO))
summary(mod3)


stargazer(mod1,mod2,mod3, type="html",out="robust3_cons.html",
          star.cutoffs=c(0.1,0.05,0.01))


